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ABSTRACT 

Context. Turbulent fluxes of angular momentum and enthalpy or heat due to rotationally affected convection play a key role in 
determining differential rotation of stars. Their dependence on latitude and depth has been determined in the past from convection 
simulations in Cartesian or spherical simulations. Here we perform a systematic comparison between the two geometries as a function 
of the rotation rate. 

Aims. Here we want to extend the earlier studies by using spherical wedges to obtain turbulent angular momentum and heat transport 
as functions of the rotation rate from stratified convection. We compare results from spherical and Cartesian models in the same 
parameter regime in order to study whether restricted geometry introduces artefacts into the results. In particular, we want to clarify 
whether the sharp equatorial profile of the horizontal Reynolds stress found in earlier Cartesian models is also reproduced in spherical 
geometry. 

Methods. We employ direct numerical simulations of turbulent convection in spherical and Cartesian geometries. In order to alleviate 
the computational cost in the spherical runs, and to reach as high spatial resolution as possible, we model only parts of the latitude 
and longitude. The rotational influence, measured by the Coriolis number or inverse Rossby number, is varied from zero to roughly 
seven, which is the regime that is likely to be realised in the solar convection zone. Cartesian simulations are performed in overlapping 
parameter regimes. 

Results. For slow rotation we find that the radial and latitudinal turbulent angular momentum fluxes are directed inward and equa- 
torward, respectively. In the rapid rotation regime the radial flux changes sign in accordance with earlier numerical results, but in 
contradiction with theory. The latitudinal flux remains mostly equatorward and develops a maximum close to the equator. In Cartesian 
simulations this peak can be explained by the strong 'banana cells'. Their effect in the spherical case does not appear to be as large. 
The latitudinal heat flux is mostly equatorward for slow rotation but changes sign for rapid rotation. Longitudinal heat flux is always 
in the retrograde direction. The rotation profiles vary from anti-solar (slow equator) for slow and intermediate rotation to solar-like 
(fast equator) for rapid rotation. The solar-like profiles are dominated by the Taylor-Proudman balance. 
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1. Introduction 

The surface of the Sun rotates differentially: the rotation pe- 
riod at the poles is roughly 35 days as opposed to 26 days at 
the equator. Furthermore, the internal rotation of the Sun has 
been revealed by helioseismology (e.g. Thompson et al. 2003): 
the radial gradient of £1 is small in the bulk of the convec- 
tion zone, whereas regions of strong radial differential rotation 
are found near the base and near the surface of the convection 
zone. According to dynamo theory, large-scale shear plays an 
important role in generating large-scale magnetic fields (e.g., 
Moffatt [19781 Krause & Radler 1980). More specifically, large- 
scale shear lowers the threshold for dynamo action and the com- 
bined effect of helical turbulence and shear yields oscillatory 
large-scale magnetic fields, resembling the observed solar ac- 
tivity pattern (e.g. Yoshimura 1 1 9751 ) . It is even possible to drive 
a large-scale dynamo in nonhelical turbulence with shear (e.g., 
Br anden burg[2S05l Yousef et al. I2008al l2008bl Brandenburg et 
al. 2008). Thus, it is of great interest to study the processes that 
generate large-scale shear in solar and stellar convection zones. 

Differential rotation of the Sun and other stars is thought 
to be maintained by rotationally influenced turbulence in their 
convection zones. In hydrodynamic mean-field theories of stel- 
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lar interiors the effects of turbulence appear in the form of tur- 
bulent fluxes of angular momentum and enthalpy or heat (cf. 
Riidiger 1989; Riidiger & Hollerbach 2004). These fluxes can be 
defined by Reynolds averaging of products of fluctuating quan- 
tities, v.i.z., the fluxes of angular momentum and heat, respec- 
tively, are 

Qij =W P (D 
Ft = cpp^T 7 . (2) 

Here overbars denote azimuthal averaging, primes denote fluc- 
tuations about the averages, Qij is the Reynolds stress, Fj is 
the turbulent convective energy flux, u is the velocity, T is the 
temperature, p is density, and cp is the specific heat at constant 
pressure. 

Much effort has been put into computing these correlations 
using analytical theories (e.g., Riidiger 1980 1982; Kitchatinov 
& Rudiger lT993l Kitchatinov et al. 119941 . Most of the analyti- 
cal studies, however, rely on approximations such as first-order 
smoothing, the applicability of which in the stellar environments 
can be contested. In order to get more insight, idealised nu- 
merical simulations, often working in Cartesian geometry, have 
been extensively used to compute the stresses for modestly large 
Reynolds numbers (e.g., Pulkkinen et al. 1993; Brummell et al. 
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1998; Chan|200TJ Kapyla et al. 2004, Rudiger et al. 2005b). 
However, the Cartesian simulations have yielded some puzzling 
results, such as the latitudinal angular momentum flux having a 
very strong maximum very close to the equator (e.g., Chan 2001 1 
Hupfer et al. 120051) and a sign change of the corresponding radial 
flux (Kapyla et al. 2004). Neither of these effects can be recov- 
ered from theoretical studies (Rudiger & Hollerbach 2004]) or 
simpler forced turbulence simulations (Kapyla & Brandenburg 
2008). The Reynolds stresses have also been computed from 
high resolution spherical convection simulations (e.g. DeRosa 
et al. 120021 Miesch et al. 2008]), but a detailed comparison with 
Cartesian results is lacking in the literature. 

Rotation also affects the turbulent convective energy trans- 
port. In fact, in the presence of rotation, the turbulent heat 
transport due to convection is no longer purely radial (e.g., 
Brandenburg et al. 119921 Kitchatinov et al. 1994; Brun & 
Rempel 2009). In a sphere, such anisotropic heat transport 
leads to latitude-dependent temperature and entropy distribu- 
tions. Such variations can be important in determining the rota- 
tion profile of the Sun: neglecting the Reynolds stress and molec- 
ular diffusion, the evolution of the azimuthal component of vor- 
ticity, ZJ = V x u, is governed by 



dt 



(3) 



where d/dz = fl ■ V is the derivative along the unit vector of 
the rotation vector, fi = (cos 8, — sin 8, 0), and p is the pressure. 
The last term on the rhs describes the baroclinic term which can 
be written as 



rep o8 



(4) 



where g = \g\is the acceleration due to gravity, s is the specific 
entropy, and VT « g/cp has been used for the adiabatic tem- 
perature gradient. In the absence of latitudinal entropy gradients, 
the solution of Eq. ([3]) is given by the Taylor-Proudman theorem, 
i.e. dfl jdz = 0. In general, however, the thermodynamics can- 
not be neglected and latitudinal gradients of entropy influence 
the rotation profile of the star via the baroclinic term. Such an 
effect is widely considered to be instrumental in breaking the 
Taylor-Proudman balance in the solar case (e.g., Rempel 2005; 
Miesch et al. 2006 ). Local simulations can be used to determine 
the latitudinal heat flux but by virtue of periodic boundaries, no 
information about the latitudinal profile of entropy can be ex- 
tracted from a single simulation. Earlier local studies suggest 
that in the presence of rotation the latitudinal heat flux is directed 
towards the poles (e.g. Rudiger et al. I2005bl ) and mean-field 
models in spherical geometry indicate that such a flux leads to 
warm poles and a cooler equator (e.g. Brandenburg et al. 1992 ), 
thus alleviating the Taylor-Proudman balance. Computing the 
turbulent heat fluxes in spherical geometry in order to compare 
with earlier results is one of the principal aims of the present 
study. Of particular importance is the sign and magnitude of the 
latitudinal heat flux. 

It is possible that the use of Cartesian geometry and pe- 
riodic boundaries give rise to artefacts which are not present 
in fully spherical geometry. In the present paper we undertake 
the computation of Reynolds stress and turbulent heat transport 
from simulations in spherical geometry as functions of rotation, 
and compare them with Cartesian simulations of the same sys- 
tem located at different latitudes. One of the most important 
goals of the paper is to find out whether the present results in 



Cartesian geometry compare with early similar studies and to 
test if these results are still valid when spherical geometry is 
used. Earlier studies comparing spherical and Cartesian mod- 
els used limited two-dimensional geometry in the spherical case 
Hupfer et. al {2006 ) whereas we perform all our simulations in 
three dimensions. Furthermore, Robinson & Chan (2001) used 
spherical wedges to compute the rotation profiles and turbulent 
fluxes using two representative runs. Here we explore a signifi- 
cantly larger portion of parameter space. As a side result we also 
obtain angular velocity profiles as a function of rotation from 
our spherical simulations which, however, are dominated by the 
Taylor-Proudman balance in the regime most relevant to the Sun. 
Thus we fail in reproducing the solar rotation profile which is a 
common problem that can currently be overcome only by in- 
troducing some additional poorly constrained terms, e.g. a lat- 
itudinal entropy gradient, by hand rather than self-consistently 
(e.g. Miesch et al. 2006). Another important use for the results 
will be the more ambitious future runs where subgrid-scale mod- 
els of the turbulent effects can be used to overcome the Taylor- 
Proudman balance. 

2. Model 

Our spherical model is similar to that used by Kapyla et al. 
(120 1 Oab but without magnetic fields. We model a segment of a 
star, i.e. a "wedge", in spherical polar coordinates where (r, 6, </>) 
denote the radius, colatitude, and longitude. The radial, lati- 
tudinal, and longitudinal extents of the computational domain 
are given by 0.65R < r < R, 8 < 8 < 180° - 8 Q , and 
< < <po, respectively, where R is the radius of the star. 
In all of our runs we take 8q ~ 15° and <fio = 90°. We study 
the dependence of the results on domain size in Appendix lAl 
In our Cartesian runs, the coordinates (x, y, z) correspond to ra- 
dius, latitude and longitude of a box located at a colatitude 8. Our 
domain spans from 0.65i? < x < R, — 0.35R < y < 0.35R and 
— 0.35i? < z < 0.35R, i.e., the extension of the horizontal di- 
rections is twice the vertical one, as has been used in previous 
Cartesian simulations (e.g. Kapyla et al. 120041 ). 

In both geometries, we solve the following equations of com- 
pressible hydrodynamics, 



D\np 
Dt 



= V ■ it, 



Du „ 1 . „ „ 

= g-2Qxu + -{V- 2vpS - Vp) 

Dt p 



Ds _ 1 
~D~t~ "pT 



( V • KVT + 2^S 2 - r cool ) 



(5) 
(6) 
(7) 



where D/Dt = d/dt + u ■ V is the advective time derivative, 
v is the kinematic viscosity, K is the radiative heat conductivity, 
and g is the gravitational acceleration given by 



g = 



GM 



(8) 



where G is the gravitational constant, M is the mass of the star, 
and r is the unit vector in the radial direction. Note that in the 
Cartesian case x corresponds to the r direction so that all radial 
profiles in spherical coordinates directly apply to the Cartesian 
model. We omit the centrifugal force in our models. This is con- 
nected with the fact that the Rayleigh number is much less than 
in the Sun, which is unavoidable and constrained by the numer- 
ical resolution available. This implies that the Mach number is 
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larger than in the Sun. Nevertheless, it is essential to have real- 
istic Coriolis numbers, i.e. the Coriolis force has to be larger by 
the same amount that the turbulent velocity is larger, but without 
significantly altering the hydrostatic balance that is determined 
by gravity and centrifugal forces. 

The fluid obeys the ideal gas law with p = (7 — l)pe, where 
7 = cp/cy = 5/3 is the ratio of specific heats in constant pres- 
sure and volume, respectively, and e = cyT is the internal en- 
ergy. The rate of strain tensor S is given by 



(9) 



where the semicolons denote covariant differentiation (see Mitra 
et al. l2009l for details). 

The computational domain is divided into three parts: a 
lower convectively stable layer at the base, convectively unsta- 
ble layer and a cooling layer at the top mimicking the effects 
of radiative losses at the stellar surface. The radial positions 
{ri,r 2 ,r 3 ,r 4 ) = (x 1 ,x 2 ,x 3 ,x 4 ) = (0.65, 0.7, 0.98, 1)R give 
the locations of the bottom of the domain, bottom and top of the 
convectively unstable layer, and the top of the domain, respec- 
tively. The last term on the rhs of Eq. (0 describes cooling in the 
surface layer given by 



r /(r) 



1 sO 



(10) 



where f(r) is a profile function equal to unity in r > r% and 
smoothly connecting to zero below, and Tq is a cooling lumi- 
nosity chosen so that the sound speed in the uppermost layer 

2 s {r = r 4 ). 



relaxes toward Cg 



2. 1. Initial and boundary conditions 

For the thermal stratification we adopt a simple setup that can be 
described analytically rather than adopting profiles from a solar 
or stellar structure model as in, e.g., Brun et al. (2004). We use a 
piecewise polytropic setup which divides the domain into three 
layers. The hydrostatic temperature gradient is given by 



dT 
dr 



-9 



cv(7 -!)(» + !)' 



(11) 



where n = n(r) is the radially varying polytropic index. This 
gives the logarithmic temperature gradient V (not to be confused 
with the operator V) as 



V = <91nT/dlnp = (n + 1) 



(12) 



The stratification is unstable if V — V a d > where V a d = 
1 — I/7, corresponding to n < 1.5. We choose n = 6 for the 
lower overshoot layer, whereas n = 1 is used in the convectively 
unstable layer. A polytropic setup with n = 1 is commonly used 
in convection studies (e.g. Hurlburt et al. l 19841 >. This implies that 
about 80 per cent of the energy is transported by radiation (cf. 
Brandenburg et al. 2005), regardless of the vigor of convection 
and the value of the Reynolds number. 

Density stratification is obtained by requiring hydrostatic 
equilibrium. The thermal conductivity is obtained by requiring 
a constant luminosity Lq throughout the domain via 



K = 



4irr 2 dT/dr 



(13) 



In order to expedite the initial transient due to thermal relaxation, 
the thermal variables have a shallower profile, corresponding to 




0.7 0.8 0.9 

r/R 



Fig. 1. Radial profiles of entropy, temperature, density, and pres- 
sure in the initial state (solid lines) and the in the saturated state 
(dashed) of Run B0. Reference values To and po are taken from 
the bottom of the convectively unstable layer in the initial state. 
The dotted vertical lines at r 2 = 0.7R and = 0.98i? denote 
the bottom and top of the convectively unstable layer, respec- 
tively. 



p oc T 1 - 4 , in the convection zone and n = 1 is only used for the 
thermal conductivity. This gives approximately the right entropy 
jump that corresponds to the required flux (cf. Brandenburg et 
al. l2.0051 >. In Fig.Q]we show the initial and final stratifications of 
specific entropy, temperature, density, and pressure for a partic- 
ular run. 

In the spherical models the radial and latitudinal boundaries 
are taken to be impenetrable and stress free, according to 



U r = 0, 



due _ ue_ du^, _ . 
dr r ' dr r 



du r _ 

w = ue = > 



-Q(f= U <t> COt ( 



ri,r 4 ), (14) 
= e 0l n-9 ). (15) 



On the latitudinal boundaries we assume that the thermodynamic 
quantities have zero first derivative, thus suppressing heat fluxes 
through the boundary. 

In Cartesian coordinates we use periodic boundary condi- 
tions in the horizontal directions (y and z), and stress free con- 
ditions in the x direction, i.e., 



du y du z 
u x = —— = —— = (x = xi,x i ). 
ox ox 



(16) 



The simulations were performed using the PENCIL Cod^3, 
which uses sixth-order explicit finite differences in space and a 
third-order accurate time stepping method (see Mitra et al. 2009 
for further information regarding the adaptation of the PENCIL 
Code to spherical coordinates). 

2.2. Nondimensional quantities 

Dimensionless quantities are obtained by setting 

R = GAI = po = c P = 1 , (17) 



http : / /pencil-code . googlecode . com/ 
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Table 1. Summary of the spherical runs. 



Run 


grid 


Ra 


Pr 




C 




Ma 


Re 


Co 


-Ether 


Ekin 




E mc r/Ekin 


E r ot/Ekin 


An/n cq 


AO 


128 x 256 x 128 


3.1 


10" 


1.0 


1.4 


10" 


4 


0.023 


38 


0.00 


0.116 


7.7 


10" 


b 


0.045 


0.004 


— 


Al 


128 x 256 x 128 


3.1 


10 b 


1.0 


1.4 


10" 


4 


0.022 


36 


0.13 


0.114 


6.9 


10" 


b 


0.016 


0.022 


-0.15 


A2 


128 x 256 x 128 


3.1 


10 6 


1.0 


1.4 


10" 


-4 


0.022 


36 


0.25 


0.114 


7.2 


10" 


-5 


0.015 


0.073 


-0.31 


A3 


128 x 256 x 128 


3.1 


10 6 


1.0 


1.4 


• 10" 




0.022 


37 


0.50 


0.113 


1.2 


10" 




0.010 


0.438 


-1.03 


A4 


128 x 256 x 128 


3.1 


10 6 


1.0 


1.4 


• 10" 


-4 


0.029 


48 


0.94 


0.112 


1.1 


10" 


_3 


0.016 


0.927 


-1.74 


A5 


128 x 256 x 128 


3.1 


10 6 


1.0 


1.4 


• 10" 


-4 


0.022 


36 


2.56 


0.111 


9.9 


10" 


4 


0.002 


0.949 


-0.37 


A6 


128 x 256 x 128 


3.1 


10 6 


1.0 


1.4 


■ 10" 


-4 


0.018 


30 


6.09 


0.114 


2.3 


10" 


4 


0.000 


0.824 


+0.20 


BO 


128 x 512 x 256 


8.6 


10" 


1.0 


1.4 


10" 


4 


0.020 


54 


0.00 


0.113 


5.8 


10" 


5 


0.036 


0.009 




Bl 


128 x 512 x 256 


8.6 


10 b 


1.0 


1.4 


■ 10" 


4 


0.020 


57 


1.34 


0.112 


6.5 


10" 


4 


0.009 


0.927 


-1.10 


B2 


128 x 512 x 256 


8.6 


10 6 


1.0 


1.4 


10" 


-4 


0.018 


50 


3.06 


0.113 


1.2 


10" 


-4 


0.001 


0.689 


+0.12 


B3 


128 x 512 x 256 


8.6 


10 6 


1.0 


1.4 


• 10" 


-4 


0.016 


44 


6.93 


0.113 


1.8 


10" 


-4 


0.000 


0.833 


+0.20 


CI 


128 x 256 x 128 


1.7 


10' 


6.7 


7.5 


10" 


5 


0.008 


12 


7.58 


0.114 


3.8 


10" 


5 


0.000 


0.817 


+0.12 


Dl 


256 x 512 x 256 


6.0 


10' 


1.0 


3.1 


• 10" 


5 


0.012 


90 


5.07 


0.113 


7.9 


10" 


5 


0.000 


0.796 


+0.20 


D2 


256 x 512 x 256 


6.0 


10 7 


1.0 


3.1 


10" 


-5 


0.012 


89 


7.68 


0.113 


1.5 


10" 


-4 


0.000 


0.895 


+0.20 



Notes. Here Ma = « rms / \/GM/R, Afl = fi cq — f2 po i c , where fl cq — tt(r4,8 = 90°) and f2 po ic = Q(r4,6 = 8 a ). Ether = (pe) and 
Ekin = (iipu 2 ) are the volume averaged thermal and total kinetic energies, respectively, in units of GMpo/R. E mer = ^(p(ug + and 
E ro t = \ (pw|) are the kinetic energies of the meridional circulation and differential rotation. 



where po is the density at r%, The units of length, velocity, den- 
sity, and entropy are then given by 



[x]=R, [u] = y/GM/R , [p] = po , [s] = c P 



(18) 



The Cartesian simulations have been arranged so that the thick- 
ness of the layers is the same, g = — (GM/x 2 )x, and R, which 
is still our unit length, has no longer the meaning of a radius. The 
simulations are governed by the Prandtl, Reynolds, Coriolis, and 
Rayleigh numbers, defined by 



Pr = 



Ra 



— , Re = —— 

Xo vk f 

GM (Ar) 4 



Co = 



2 fin 



vxoR 2 



1 ds 
cp dr 



(19) 
(20) 



where xo = K/(p m cp) is the thermal diffusivity, fef = 2n/Ar 
is an estimate of the wavenumber of the energy-carrying eddies, 
Ar = r3 — r2 is the thickness of the unstable layer, p m is the 
density in the middle of the unstable layer at r m = (7-3 + 

and u rms = ^J^{u 2 . + uf) is the rms velocity, where the angular 

brackets denote volume averaging. In our definition of it rms we 
omit the contribution from the (^-component of velocity, because 
it is dominated by the large-scale differential rotation that devel- 
ops when rotation is included. The entropy gradient, measured 
at r m in the initial non-convecting state, is given by 



1 ds\ _ V,_ 
cp dr ) r Hp 



V 



ad 



(21) 



where V m = (d\nT/dlnp) rra , and Hp is the pressure scale 
height at r m . 

The energy that is deposited into the domain at the base is 
controlled by the luminosity parameter 



c = 



Po {GMf/ 2 Rin 



(22) 



where Lq = 4irrlF\> is the constant luminosity, and — 
— (KdT/dr)\ r=ri is the energy flux imposed at the lower 



\ 
-5 



1.0 



0.5 



0.0 



■ 1 1 1 I 1 1 * 1 1 1 1 1 1 1 1 1 1 1 1 1 ' 1 1 1 1 1 1 1 1 1 1 ! 1 
■ / il.. 1 ■ 




' T7 

'■■ / \ " 




V : 




A i 


/ 

: 1 
: 1 


/ \ 

"~ - ^ ^cool / 




i .... 1 .... 1 .... 1 . . 





0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 
r/R 

Fig. 2. Radiative (dotted line), enthalpy (dashed), kinetic energy 
(dash-dotted), cooling (long dashed), and viscous (triple-dot- 
dashed) luminosities as functions of radius from Run AO. The 
solid line shows the sum of all fluxes, and the red dashed line the 
luminosity Lq fed into the domain through the lower boundary. 
The vertical dotted lines at r = 0.7R and r = 0.98i? denote re- 
spectively the bottom and top of the convectively unstable layer 
in the initial state. 

boundary. Furthermore, the stratification is determined by the 
pressure scale height at the surface 



(7 - i)c v r 4 

GM/R 



(23) 



where T4 = T(r = rjj. Similar parameter definitions were used 
by Dobler et al. (120061 1 . We use £ = 0.020, which results in a 
density contrast of 10 2 across the domain. 

3. Results 

Our main goal is to extract the turbulent fluxes of angular mo- 
mentum and heat as functions of rotation from our simulations. 
In order to achieve this we use a moderately turbulent model and 
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Fig. 3. Radial velocity u r at a small distance (r = 0.9R) below the surface from Runs B0-B3. The scales give u r in units 
of the local sound speed. For visualization purposes, the domain is duplicated fourfold in the longitudinal direction. See also 

http : / / www .helsinki . fi/ ~kapy la /movies . html 



vary the rotation rate, quantified by the Coriolis number, from 
zero to roughly six in Set A (see Table [TJ. We also perform a 
subset of these simulations at higher resolution in Set B and a 
three runs (CI, Dl, and D2) with a lower Mach number. The 
runs in Set A were initialized from scratch, whereas in Set B 
a nonrotating simulation BO was run until it was thermally re- 
laxed. The runs with rotation (B1-B3) were then started from 
this snapshot and computations carried out until a new saturated 
state was reached. The runs Dl and D2 were remeshed from a 
non-rotating, thermally relaxed model at a lower resolution. In 
Fig. Q] we compare the initial and final stratification of specific 
entropy, temperature, density, and pressure for Run BO. 

As noted in Sect. 12.11 our polytropic setup leads to a sys- 
tem where radiative diffusion transports 80 per cent of the total 
energy. We show the flux balance in the statistically saturated 
state from Run AO in Fig. [2] where the different contributions 
are given in terms of luminosities Li = Airr 2 Fi, and where 



completely defined by I, m, and v, ™ since, assuming approxi 





dr ' 


(24) 


± conv — 


-c P pu' r T', 


(25) 




^pu 2 u r , 


(26) 


F ■ — 

± vise — 


-2upUiSi r , 


(27) 


-^cool — 


J r coo i dr. 


(28) 



Here we consider averages over (ft and 9. We find that in the non- 
rotating case the convective flux accounts for roughly 30 per cent 
of the total luminosity and the (inward) kinetic energy flux is be- 
tween 10 and 15 per cent. When rotation is increased, both the 
convective and kinetic fluxes decrease. The viscous flux is al- 
ways negligible. The cooling flux transports the total luminosity 
near the surface. 

Visualizations of u r at a small distance below the surface are 
shown in Fig. [3] for Runs B0-B3. The convective velocities u' 
can be decomposed in terms of poloidal (u' P ) and toroidal (u' T ) 
parts following Lavely & Ritzwoller ( 1 19921 1 



v!p = Real V {4 m (r)r™f + 4 m (r)Vl^ m } 



I. I'll 



(30) 



where Y [ m (8, (ft) are spherical harmonics of degree I and order 
rn. The geometry and amplitude of the poloidal velocity are 



mate mass conservation, uK" and wff 1 are related as 



,,lrn 



(r) 



d r (r 2 pu l ™(r)) 
prl{l + 1) 



(31) 



The poloidal flow has characteristics of Benard convection cells 
with upwellings at the centres of cells and downdraughts on the 
peripheries. The toroidal flows are characterised by their ampli- 
tude and geometry given by wlf 1 , I, and m respectively. In con- 
trast to poloidal flows, their nature resembles that of rotation, 
jets or horizontal vortices. In Fig. [3] we observe that so called 
banana cells become prominent in the radial velocity with an in- 
crease in the Coriolis number. Such structures are poloidal flows 
given by spherical harmonic Y^ n {9 1 (ft). For Run B3 in Fig. [3] 
we find maximum power at m = 16. Note that the reality of the 
banana cells in the Sun is hotly debated. Even though significant 
power is found at wavenumbers corresponding to giant cells in 
the surface velocity spectra of the Sun, no distinct peak has been 
found at those wavenumbers (Chou et al. |199U Hathaway et al. 
2000). Global helioseismology caps the maximum radial veloc- 
ity of the banana cells at 50 m s _1 (Chatterjee & Antia 2009). 
We study the importance of the banana cells to the Reynolds 
stresses in more detail in Sects. l3~Tl and l3.2.T1 

3. 1. Reynolds stress 

The angular momentum balance of a star is governed by the con- 
servation law (Rudiger ll989l l 



0_ 

Of 



{pm 2 n) = -v 



pvj I VJ 



u'.u- 



(32) 



where m = r sin 9 is the lever arm and Tt mcr = (u r ,ug) is the 
meridional circulation. The latter term on the rhs describes the 
effects of the Reynolds stress components and Qg<p, which 
describe radial and latitudinal fluxes of angular momentum, re- 
spectively. The stress is often parameterised by turbulent trans- 
port coefficients that couple small-scale correlations with large- 
scale quantities, i.e. 



(29) Q y = Ayfcfifc - Mi 



du k 

OXi 



(33) 



where A^fe describes the nondiffusive contribution (A-effect) 
and Mijki the diffusive part (turbulent viscosity), cf. Riidiger 
( 1989). However, disentangling the two contributions is not pos- 
sible, see e.g., Snellman et al. d2Q09b and Kapyla et al. (2010b). 
We postpone a detailed study of the turbulent transport coeffi- 
cients to a future study and concentrate on comparing the total 
stress with simulations in Cartesian geometry. 
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Fig. 4. Vertical Reynolds stress, Q r ^, from Set A. 

It is convenient to display the components of the Reynolds 
stress in non-dimensional form (indicated by a tilde), and to de- 
fine 

Qij =Ui u j/ u tm S , (34) 

where u lms — u lms (r,9) is the meridional rms-velocity. The 
averages are calculated over the azimuthal direction and time 
also for u rms . In the following, we refer to the three off-diagonal 
components, Q r ^, Qg^, and Q r g, as vertical, horizontal, and 
meridional components, respectively. Representative results for 
the vertical stress component Q r< ^ are shown in Fig. |4] We find 
that for slow rotation (Run Al), Q rt p is small and does not ap- 
pear to show a clear trend in latitude. In Run A2 with Co w 0.25 
the stress is more consistently negative within the convectively 
unstable layer, showing a symmetric profile with respect to the 
equator. These two runs tend to show the largest signal near the 
latitudinal boundaries which is most likely due to the boundary 
conditions there. Similar distortions are also seen in the large- 
scale flows (see Sect. 13.41 ). In the intermediate rotation regime 
(Runs A3-A5), Q r ^ is predominantly negative, although regions 
of opposite sign start to appear near the equator. In Run A6 
the stress is mostly positive. Qualitatively similar results are ob- 
tained from the runs in Set B, Runs CI, Dl, and D2. Therefore 
there is a sign change roughly at Co = 2. The results for most 
quantities from Runs B2 and Dl with intermediate values of 
Co are similar to those of Runs A5 and A6, respectively. Thus, 
we usually show results only from Runs A4, A5, and A6 in or- 
der to demonstrate the qualitative change that occurs for many 
quantities in the range Co w 1 ... 6. A similar phenomenon has 
been observed in Cartesian simulations (Kapyla et al. 2004}. We 
note that the behaviour of Q r( f, in the most rapidly rotating runs, 
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Fig. 5. Horizontal Reynolds stress, Qg,p, from Set A. 

namely a small negative region at the equator and a positive peak 
near the surface at somewhat higher latitudes was also reported 
by Robinson & Chan (EOOTi 

We find that the horizontal stress, Qg^, is almost always 
positive (negative) in the northern (southern) hemisphere for 
Co < 1, i.e. antisymmetric about the equator, see Fig. [5] For 
intermediate rotation (Runs A4 and A5) the stress is observed 
to change sign at high latitudes. In Fig.|6]we plot the latitudinal 
profiles of the horizontal stress and the mean angular velocity at 
different depths for the Runs A4-A6. It can be seen that near the 
bottom of the convection zone, the profile of the stress becomes 
more and more concentrated about the equator as the Coriolis 
number increases. An especially abrupt change can be observed 
for Run A5 (Co w 2). A similar peak also persists in Runs A6, 
B3, CI, and D2 with the largest Coriolis numbers. Note, how- 
ever, that the sign of the latitudinal differential rotation changes 
as Co increases to six for Run A6. The results of Robinson & 
Chan (2001) also show a peak of Qg^, occurring at a latitude 
range 10° . . . 15°, depending on depth. 

Using Eqs. d29b-d30l). we can calculate the stress Qg^ = 
J2i,y,m Qe ™ b y azimuthal averaging, with 

rfl'va _ 1 Am ...I'm ( 1 dPF dPjP r jf_ pmpm \ nS x 

where P™ (0) are the associated Legendre polynomials and xn = 
r sin 0. Note that I and I' denote the degrees of the poloidal and 
the toroidal flow, respectively. It is easy to see that the contribu- 
tion to the azimuthally averaged Qg^ is always zero from cross- 
correlation between two poloidal velocity fields. Finite contribu- 
tions to Qg^ instead come from correlations between poloidal 
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Fig. 6. Latitudinal profiles of Qg^ and Q for Runs A4, A5, and A6 (from left to right) at three different depths (solid 0.7 R, dotted 
0.8i?, dashed 0.9R). The red dotted lines show data from corresponding Runs Bl, B2, and B3 from r = 0.8R. The solid red lines 
in the right panels show data from Run D2 at r = 0.8R. The open red diamonds in the top panels denote Cartesian Runs cAl- 
cA4, cDl-cD4, and cEl-cE4, from left to right. The blue squares in the top-middle panel show the values of Q yz computed from 
Fourier-filtered velocity fields from Runs cDl-cD4. Note that only a part of the full latitudinal range is shown. 



flow and toroidal flow having the same azimuthal degree m. We 
have used small-scale velocity fluctuations (i.e., m ^ modes) 
to calculate the Reynolds stresses in the numerical simulations 
according to Eq. d34l i. The finite correlation of the rotation and 
the meridional flow are not included in this discussion since both 
are characterised by to = and thus do not correspond to our 
definition of velocity fluctuations. 

Recently, Bessolaz & Brun (1201 11 ) have used wavelets and 
autocorrelation techniques to unravel the structure of giant cells 
in their 3 -dimensional hydrodynamic convection simulations. It 
is an involved exercise to calculate the net stress by estimating 
the power in each triplet (/, /', m) by wavelet analysis. It is, how- 
ever, possible to look for certain combinations of Legendre poly- 
nomials that can contribute to the peaks of Qg^ near the equa- 
tor as obtained from numerical simulations in spherical geome- 
try. A visual inspection of the radial flows in Fig.[3]for Run B2 
shows four prominent banana cells within the domain which ex- 
tends from to 7r/4 in the azimuthal direction, which means 
that the angular dependence is most likely Y{q , Hence we set 
/ = 16, 1' = 16, 17 for the calculation of the stresses and vary 
to in search for a match between the peaks of Qg^ from the runs 
A1-A6 and Eq. (l3Sb . We illustrate the angular part of Qg'^ 11 , 
for particular values of I, I' and to in Fig. [7] We can see from 
here that peaks in Q 1 ^ 17:15 (dashed line) appear at ±6° as well 
as at ±20° latitude, whereas peaks in Qg^' 17 ' 16 appear at ±10° 
latitude, and the highest peaks in Q 1 ^ 17:8 appear at ±60° lati- 
tude. Comparing Fig. [5] with Fig. O we see that at slow rotation 
(Runs Al and A2), a major contribution to the stress may come 
from giant cells with an angular dependence Y^. At higher Co, 
the stress may have contributions from banana cells with angular 
dependence (compare solid line in top right panel of Fig. [6] 



with solid line of Fig. [7]). We shall return to the question regard- 
ing the contribution of banana cells in the context of Cartesian 
runs in Sect. 13.2.11 However there also exists symmetric con- 
tribution to Qg,p from components like Qgl 16 ' 16 , but we do not 
see any significant symmetric part in the horizontal stresses from 
the numerical simulations. On this basis, zonal flows of the form 
wi^r x Vij' can be said to be negligible in spherical convection 
simulations. These zonal flows correspond to a row of horizontal 
vortices with their centres on the equator. 



Finally, let us discuss the stress component Q r g. It does 
not directly contribute to angular momentum transport, but it 
can be important in generating or modifying meridional circu- 
lation, and it has routinely been considered also in earlier stud- 
ies (e^, Pulkkinen et al. 119931 Rieutord et al. 119941 Kapyla et 
al. 2004). Figure 8 shows the stress component Q r g from Set A. 
We find that for slow rotation (Run Al) the stress is quite weak 
and shows several sign changes as a function of latitude. It is 
not clear whether this pattern is real or an artefact of insuffi- 
cient statistics. For intermediate rotation (Runs A2-A4), Q r g 
shows an antisymmetric profile with respect to the equator being 
positive in the northern hemisphere and negative in the south, 
in accordance with earlier Cartesian results (e.g. Kapyla et al. 
2004). Although the theory for this stress component is not as 
well developed as that of the other two off-diagonal components, 
Riidiger et al. (2005a) state that Q r g should always be negative 
in the northern hemisphere, which is at odds with our results. 
However, in our rapid rotation models (Runs A5-A6) the sign is 
found to change. 
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Fig. 7. Angular part of Q 1 ^™ normalized by the maximum value 
for four different cases characterized by triplets (I, I', m) as in- 
dicated by the legend. The latitudes of the peaks for the triplets 
are indicated on the respective curves. 

3.2. Comparison with Cartesian simulations 

Before describing the Reynolds stress obtained from our sim- 
ulations in Cartesian coordinates, we note that the rms veloci- 
ties in the Cartesian runs are in general almost twice as large as 
in the spherical ones with the same input parameters (compare, 
e.g., Run AO in Table Q] and Run cAO in Table [2]). We argue in 
Sect. l3.3l that this is the result of adopting a radial dependence of 
gravity in the plane-parallel atmosphere. 

The radial profiles of the three off-diagonal components of 
the Reynolds stress in Cartesian coordinates agree with previous 
studies (Kapyla et al. 120041 Hupfer et al. 120051) for the range 
of latitudes and Coriolis number explored here (compare Fig. [9] 
with bottom panel of Fig. 11 of Kapyla et al.[2004 and Figs. 3 
and 5 of Hupfer et al. I20051 >. For moderate rotation (Runs cAl- 
cA4), the vertical component Q xz (left panels of Fig. [9} is neg- 
ative in the bottom part of the convection zone and almost zero 
at the top. The cases with Co 2.3 (Runs cDl-cD4) show neg- 
ative values at the bottom and positive values at the top of the 
convection zone. For Co « 4.0 (Runs cEl-cE4), the amplitude 
of the positive part of the stress near the surface increases and 
the negative part at the bottom decreases. We notice that the spa- 
tial distribution of Q xz , as well as its variation with the Coriolis 
number, are in a fair agreement with the corresponding spheri- 
cal runs in the same range of Co (Runs A3-A5). In the spher- 
ical Run A6 with the highest Coriolis number of roughly six, 
the stress is observed to become predominantly positive in the 
convection zone. This is not seen in the Cartesian counterparts 
that reach Coriolis numbers of roughly four (Runs cEl-cE4), in 
which the negative peak near the bottom still persists, although 
it has decreased in magnitude. The difference is possibly due to 
the lower Coriolis number in the Cartesian runs. It is noteworthy 
that also the symmetry of this stress component with respect to 
the equator is captured by the Cartesian simulations. 

Radial profiles of the horizontal stress, Q yz , from the 
Cartesian simulations are shown in the middle panels of Fig. [9] 
and latitudinal profiles in Fig.|6]with open squares and diamonds. 
Similarly as in the spherical runs, this component peaks both at 
top and bottom of the convective layer. However, some discrep- 
ancies are observed between the profiles in different geometries. 
For instance, in spherical Run A4 the stress is somewhat more 
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Fig. 8. Meridional Reynolds stress, Q r o, from Set A. 



widely distributed than in the corresponding Cartesian runs. In 
spherical Run A5 the radial profile of this component exhibits a 
bump at the bottom of the convection zone which is much larger 
than in the corresponding Cartesian cases. Note, however, that in 
Fig. [9] the uppermost peak moves inwards with increasing rota- 
tion between Sets cA and cD, and at the same time as the lower- 
most peak increases in amplitude. For the spherical Run A6 with 
the highest Coriolis number of roughly six, the stress changes 
sign in the region near the surface, which is not visible in the 
Cartesian simulations with Coriolis numbers of roughly four 
(Runs cEl-cE4). 

Finally, the meridional Reynolds stress, Q xy , corresponding 
to Q r e, is positive in the entire convection zone for moderate 
rotation (Runs cAl-cA4). For larger Co, Q xy is negative in the 
lower part of the domain (see the right panels of Fig. |9). Similar 
behaviour occurs in the spherical case with intermediate rotation 
(Runs A3-A5). In the most rapidly rotating case (Run A6) an- 
other sign change occurs near the equator (see Fig.O, which is 
not observed in Cartesian runs. This, however, could again be 
explained by the smaller Co in the Cartesian runs. 

3.2.1. Filtering banana cells 

The large amplitude of the horizontal Reynolds stress, peaking 
around ±7° latitude, has been an intriguing issue for several 
years (e.g., Chan |200"T1 Hupfer et al. I2U531 120061 >. One factor 
that might be contributing to the Reynolds stress are the large- 
scale banana cell-like flows that develop near the equator (e.g., 
Kapyla et al. 2004; Chan 2007 ). Such flows vary in the azimuthal 
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Fig. 9. From left to right: radial profiles of Q xz , Q yz , and Q xy from Cartesian Runs cAl-cA4 (top panels), Runs cDl-cD4 (middle 
panels), and Runs cEl-cE4 (bottom panels). The red diamonds correspond to the radial profiles of the stresses in the spherical 
Runs A4-A6. The blue squares in the middle panel show Fourier- filtered data from Run cD2. 



(z) direction and can lead to overestimation of the contribution 
of turbulence, especially if averaging is performed over the az- 
imuthal (2) direction. We explore this possibility by filtering out 
the contribution coming from the large-scale structures observed 
in the yz-plane (the so-called banana cells observed in spherical 
simulations). The procedure used in this analysis is described 
below. 

We perform a Fourier decomposition of the horizontal ve- 
locities and find out at which Fourier mode the contribution of 
the large scales peaks in the spectra. We find that the maximum 
is usually situated at wavenumber q = 2. Next we remove this 
mode from the spectra and make an inverse Fourier transforma- 
tion, thus obtaining the velocity field without the contribution 
from the large-scale motions. Finally, we compute Q yz from the 
filtered velocities. 

Horizontal stress Q yz computed from filtered velocity fields 
for Runs cDl-cD4 for different latitudes at r = 0.9R are plot- 
ted with blue square symbols in Fig. [6] The radial variation of 
Q yz at 7° for Run cD2 is shown with blue square symbols in 



Fig. [9] It is clear from these figures that the q = 2 mode is the 
dominant contribution to Q yz near the surface and it also affects 
significantly the secondary peak in deeper layers. Thus, a flatter 
profile in latitude with a reduced amplitude of the stress is ob- 
tained in comparison to the non-filtered values. The maximum, 
however, still resides around ±7°, which is at odds with theory 
(e.g. Riidiger & Kitchatinov l20071 l. 

3.3. Turbulent heat transport 

In non-rotating convection the radial heat flux, 



F r = c P pu' r T', (36) 

transports all of the energy through the convection zone. 
According to mixing length theory, velocity and temperature 
fluctuations are related via u'^ ~ (AT /T)g£, where £ is the mix- 
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Table 2. Summary of the runs in Cartesian coordinates. 



Run Latitude 


Re 


Co 


Ma 


E k 






E T ot 1 -Ek 


cAO 


- 


63 


0.00 


0.038 


1.7 


■ 10" 


4 


0.052 


0.001 


cFO 


- 


28 


0.00 


0.027 


2.3 


■ 10" 


-4 


0.021 


0.001 


cFl 


- 


11 


0.00 


0.022 


2.5 


• 10" 


-4 


0.002 


0.000 


cAl 


0° 


64 


0.85 


0.039 


2.9 


• 10" 


4 


0.001 


0.288 


cA2 


7° 


65 


0.84 


0.039 


2.0 


■ 10" 


-4 


0.021 


0.017 


cA3 


14° 


65 


0.84 


0.039 


1.8 


• 10" 


-4 


0.014 


0.007 


cA4 


21° 


65 


0.85 


0.039 


1.8 


■ 10" 


4 


0.012 


0.008 


cBl 


0° 


61 


1.49 


0.037 


4.5 


• 10" 


4 


0.000 


0.623 


cB2 


7° 


70 


1.30 


0.042 


2.4 


■ 10" 


4 


0.023 


0.012 


cB3 


14° 


68 


1.33 


0.041 


2.0 


■ 10" 


-4 


0.012 


0.007 


cB4 


21° 


68 


1.34 


0.041 


1.9 


• 10" 


-4 


0.005 


0.009 


cCl 


0° 


60 


2.14 


0.036 


2.8 


• 10" 


4 


0.000 


0.347 


cC2 


7° 


76 


1.68 


0.046 


2.5 


■ 10" 


-4 


0.029 


0.031 


cC3 


14° 


72 


1.77 


0.044 


2.2 


■ 10" 


4 


0.013 


0.011 


cC4 


21° 


72 


1.78 


0.043 


2.1 


■ 10" 


-4 


0.004 


0.011 


cDl 


0° 


69 


2.38 


0.042 


7.5 


■ 10" 


4 


0.000 


0.584 


cD2 


7° 


78 


2.09 


0.047 


2.5 


■ 10" 


-4 


0.029 


0.018 


cD3 


14° 


47 


2.32 


0.043 


2.0 


■ 10" 


-4 


0.009 


0.013 


cD4 


21° 


70 


2.36 


0.042 


2.1 


• 10" 


-4 


0.003 


0.005 


cEl 


0° 


50 


3.66 


0.045 


1.2 


• 10" 


3 


0.000 


0.685 


cE2 


7° 


36 


4.00 


0.041 


1.6 


■ 10" 


-4 


0.025 


0.009 


cE3 


14° 


34 


4.24 


0.039 


1.5 


■ 10" 


-4 


0.005 


0.005 


cE4 


21° 


31 


4.67 


0.035 


1.3 


■ 10" 


-4 


0.001 


0.008 



Notes. Here, we use a resolution of 64 x 128 grid points. For the sets of 
Runs cA-cD, Ra « 3. MO 6 , and for the set of Runs cE, Ra « 1.4-10 6 . 
Thermal energy in all of the cases is Ether ~ 0.117. All quantities are 
computed using the same definitions and normalization factors as in 
Table □ 
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Fig. 10. Normalized radial turbulent heat flux raised to the 2/3 
power as a function of r (x) (solid lines). The dashed and dot- 
dashed lines correspond to the squares of the radial velocity and 
temperature fluctuations scaled with the coefficients k u and kr, 
respectively. The upper (red), middle (blue) and lower (black) 
curves correspond to Runs cAO, cFO and AO, respectively. 



ing length, gl — c 2 , and AT = V T' 2 . Thus, the three quantities 
are related via: 



AT 



,'■2 



El 

p4 



2/3 



(37) 
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Fig. 11. Turbulent heat conductivity xt from Runs AO (solid 
line) and B0 (dashed line). The inset shows the radial heat flux 
F r (solid line) and an analytical expression given in Eq. (HTI) 
(dashed line) normalized by the heat flux at r\ from Run AO. 

These quantities are shown in Fig. [I0]for non-rotating simula- 
tions in Cartesian (Run cAO) and spherical (Run AO) geometries. 
Here we use the coefficients 



cz 



« 2 /c s 2 ) 



kn 



CZ 



(AT/T) 



(38) 



where (.)cz denotes an average over the convection zone. For 
both geometries we obtain k u k, 0.4 and kx ~ 1.3, values that 
are in good agreement with previous results (Brandenburg et al. 
2005 ). Note, however, that the magnitude of the flux in Cartesian 
coordinates is around four times larger than that in the spherical 
one, implying a difference of 4 1 / 3 « 1.6 in the radial velocities 
according to Eq. J38i . This is roughly the same factor seen in 
the rms velocities (compare Runs AO and cAO). This difference 
arises from the fact that we are considering a depth dependent 
gravity also in the Cartesian simulations. In spherical geometry, 
the luminosity is constant and the flux decreases outwards pro- 
portional to r~ 2 , whereas in Cartesian geometry the flux is con- 
stant. This means that for the same profile of thermal conductiv- 
ity, a significantly larger portion of the energy is transported by 
convection in the Cartesian case. We verify this result with a sep- 
arate Cartesian model in which the radiative flux is constant and, 
like in the other models, the gravity varies with depth. In this 
case the thermal conductivity varies with radius. The profiles of 
the quantities depicted in Fig.[10]obtained from this run (see blue 
lines and Run cFO in Table |2]i are in better agreement with the 
spherical case. Similar results have been obtained if both, radia- 
tive flux and gravity, are constant (Run cFl). 

The radial turbulent heat transport may also be described in 
terms of a turbulent heat conductivity (e.g. Riidiger 1989) 



F r = cppu' r T' = -p7xtV r s, (39) 
from which we can solve the turbulent heat conductivity as 



Xt 



cpu' r T' 
TV r s ' 



(40) 



The result, normalized by a reference value xto = u rms / (3fcf ), 
for Runs AO and B0 are shown in Fig. [TT] Here averages over 
longitude and latitude are considered. We find that the value of 
Xt is almost ten times the reference value. The apparently large 
value is most likely due to the normalization factor which is 
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Fig. 12. Off-diagonal component \9r of the turbulent heat con- 
ductivity according to Eq. ( l44b from Runs Al (solid line), 
A3 (dashed), A6 (dot-dashed), B3 (triple-dot-dashed), CI (red 
dashed), and Dl (blue dotted). 



based on a volume average of the rms velocity and a more or less 
arbitrary length scale fcjT 1 (see also Kapyla et al. 1201 0b). The 
sharp peaks and negative values of \t towards the bottom and 
top of the convectively unstable region reflect the sign change of 
the entropy gradient which is not captured by Eq. (T40b . 

According to first-order smoothing (e.g. Riidiger 1989]), the 
radial flux can be written as 



F (FOSA) = 



r c ul p TV r s, 



(41) 



where r c is the correlation time of turbulence. We compare the 
actual radial heat flux with the rhs of Eq. (HTl) in the inset of 
Fig. [TT| where r c is used as a fit parameter. A reasonable fit 
within the convection zone is obtained if the Strouhal number 



St — T c u [ms k{ , 



(42) 



is around 1 .6 which is consistent with previous results from con- 
vection (e.g. Kapyla et al. 2010b). Note that the ratio Xt/xw 
gives a measure of the Strouhal number because in the general 
case Xto = ^ T cU^ ms = StM rms /(3/cf), whereas in the main panel 
of Fig. (fTTb we assume St = 1. 

In rotating convection, Eq. ( l39t no longer holds and the heat 
flux becomes latitude-dependent. In mean-field theory this can 
be represented in terms of an anisotropic turbulent heat conduc- 
tivity (Kitchatinov et al. 1994) 



Xij = XtSij + Xn£ijk^k + XnnQift. 



j i 



(43) 



where and are the Kronecker and Levi-Civita tensors 

and Cli is the unit vector along the zth component of f2. This 
indicates that non-zero latitudinal and azimuthal heat fluxes are 
also present in rotating convection. However, in order to com- 
pute all relevant coefficients from Eq. ( f43b . a procedure similar 
to the test scalar method (Brandenburg et al. 2009 ) would be re- 
quired in spherical coordinates. In most of our runs, however, 
the radial gradient of entropy is greater than the latitudinal one. 
Thus we can approximate the latitudinal heat flux by 



Fg = -pTxBr^rS - pTxee^gs w -pTxer^rS, 



(44) 



from which the off-diagonal component xer can be computed 
in analogy to Eq. d40) . Note that the sign of X9r gives the di- 
rection of the latitudinal heat flux so that positive (negative) val- 
ues indicate equatorward (poleward) in the northern (southern) 
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Fig. 13. Turbulent heat fluxes F r (top panel), Fg (middle), and 
F^ (bottom) from the runs indicated in the legend in the top 
panel. The symbols included in the top and middle panels cor- 
respond to vertical and latitudinal fluxes from Runs cAl-cA4 
(blue squares) and cEl-cE4 (red diamonds) scaled down by a 
factor of four (see the text for details). The data for Runs CI and 
D2 are scaled up by a factor of four. 



hemisphere. According to Eqs. d43l and d44l >. Fg oc sin 9 cos 9, 
indicating a sign change at the equator. 

Representative results from Runs Al, A3, A6, B3, CI, and 
Dl are shown in Fig. 112) . For slow rotation (Run Al), xer is 
small and shows no coherent latitude dependence. In the inter- 
mediate rotation regime (Run A3), xsr is positive (negative) in 
the northern (southern) hemisphere. In the most rapidly rotat- 
ing case (Runs A6 and B3), the sign changes so that the heat 
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Fig. 14. Specific entropy in the upper part of the convectively unstable layer in Runs Bl (left), Dl (middle), and D2 (right). The 
0-extent is duplicated fourfold for visualization purposes. 
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Fig. 15. Top row: radial profiles of entropy from six colatitudes as indicated by the legend in the leftmost panel from Runs Al (left 
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row at three radial positions indicated by the legend in the left panel. The red and blue dashed curves in the lower right panel show 
data at r — 0.9R from Runs CI and D2, respectively. 



flux is towards the poles. Qualitatively similar results are ob- 
tained from rapidly rotating Runs CI and Dl with a lower Mach 
number. The smoother latitude profile of \re in Run CI reflects 
the smoother entropy profile (see Fig. [T3} , The qualitative be- 
haviour as a function of rotation is similar to that found in local 
simulations (Kapyla et al. 12004). Comparing with Fig. fTTI we 
find xer/xt = Xdr/Xrr ~ 0.1, which is of the same order of 
magnitude as in local convection models Kapyla et al. (2004) 
and forced turbulence Brandenburg et al. (2009). We note that 
the latitudinal entropy gradient, which we neglected in Eq. d44l) . 
can become comparable with the radial one in the rapid rotation 
regime near the equator. Since V#s < in the northern hemi- 
sphere (cf. Fig. [15}, the latter term in Eq. ([44} yields a positive 
contribution to the flux. Thus our values of X6r near the equator 
are likely to be underestimated in the rapid rotation regime. We 
postpone a more detailed study of the turbulent transport coef- 
ficients to a future publication and discuss the different compo- 



nents of the turbulent heat fluxes. We present the components of 
convective energy flux as 

Ft = Fjpc; 3 , (45) 

where longitudinal averages are used. 

Figure Q~3] shows the normalized turbulent heat fluxes as 
functions of latitude from five runs with slow (Run Al), in- 
termediate (Run A4), and rapid (Runs A6, CI, and D2) rota- 
tion. We find that F r shows little latitudinal variation except 
near the latitudinal boundaries for slow and moderate rotation 
(Runs A1-A3). For intermediate rotation F r peaks at mid lati- 
tudes (Runs A4-A5) whereas in the most rapidly rotating cases 
(Runs A6, CI, and D2) the maxima occur near the equator and 
at the latitudinal boundaries. This behaviour follows the trend 
seen in the entropy profile (Fig[T5}: the radial gradient of entropy 
shows only a minor variation as a function of latitude in the most 
slowly rotating runs (A1-A3). In Runs A4 and A5 the gradient 
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is the steepest at mid latitudes and at the equator in Run A6. 
We find that the entropy gradient can become positive at certain 
latitudes, e.g. close to the pole for Run A4 and around latitudes 
±30° in Run A6. 

The horizontal fluxes, Fg and are negligibly small in 
comparison to the radial flux F r in the slow rotation regime 
(Run Al). The latitudinal flux is consistent with zero for all 
depths in Run Al (see Fig. [T3T >. For intermediate rotation 
(Runs A2-A4) the latitudinal flux is mostly equatorward. For the 
most rapidly rotating cases the sign changes so that in Runs A6, 
CI, and D2, F y is mostly poleward in the convection zone. The 
magnitude of the latitudinal flux also increases so that the max- 
imum values, which are located near the surface, can become 
comparable with the radial flux. The azimuthal flux is also small 
and always negative, i.e., in the retrograde or westward direc- 
tion, in accordance with the results of Rudiger et al. (20053 and 
Brandenburg et al. ([2009b . 

In some of the panels in Fig. [13] we also present results from 
Cartesian simulations (see the red and blue symbols) from the 
same depth. As discussed above, the fluxes are larger in this ge- 
ometry, due to which we have scaled the fluxes down by a fac- 
tor of four in this figure. We find that the latitude profiles of 
the radial and latitudinal heat fluxes in the Cartesian simulations 
are in rather good agreement with the spherical results. This is 
more clear in the rapidly rotating cases cEl-cE4 in comparison 
to Run A6 (see the right panels of Fig. [Oi l, where the large peak 
of F r at the equator, and the sharp peak of Fg at low latitudes 
are reproduced. 

We find that the latitudinal entropy profiles show a local 
maximum (slow and intermediate rotation) or a minimum (rapid 
rotation) at the equator, see Fig. [14] and the bottom panels of 
Fig. [15] The entropy profiles in the most rapidly rotating simu- 
lations (Run A6 and B3) are similar to that obtained by Miesch 
et al. (2000) but differs from the more monotonic profiles of e.g. 
Brun et al. (2002) and the lower Mach number case Run CI. 

3.4. Large-scale flows 

The rotation profiles from the runs in Set A are shown in Fig. [16] 
For slow rotation (Runs A1-A2), a clear large-scale radial shear, 
almost independent of latitude, develops. This is an old result 
going back to Kippenhahn ( 1963) that is expected for turbulence 
whose vertical motions dominate over horizontal ones (Rudiger 
1989). Such a result has been obtained in many mean-field mod- 
els (e.g. Brandenburg et al. 1990) and simulations since then 
(Brun & Palacios 2009 who refer to such flows as shellular). 
However, the f2-profiles in these runs are clearly different at 
high latitudes, which is probably an artefact due to the latitu- 
dinal boundaries. As the Coriolis number is increased, the ra- 
dial shear remains negative, equatorial deceleration grows, and 
the isocontours of ft tend to align more with the rotation vector 
(Runs A3-A4) - in accordance with the Taylor-Proudman theo- 
rem. Similar anti-solar rotation profiles have been reported also 
by Rieutor d et al . ([19941 Dobler et al. (I2006I Brown d2009l 
and Chan (2010). Such rotation profiles are usually the result 
of strong meridional circulation (Kitchatinov & Rudiger 2004) 
which is consistent with the present results. Run A5 represents 
a transitory case where bands of faster and slower rotation ap- 
pear, whereas in Run A6 a solar-like equatorial acceleration is 
seen. Similar transitory profiles have recently been reported by 
Chan (2010). The rotation profile in Run A6 is dominated by the 
Taylor-Proudman balance and the latitudinal shear is concen- 
trated in a latitude strip of ±30° about the equator. Similar Sl- 
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Fig. 16. Azimuthally averaged flows from the runs in Set A. The 
contours show f2 = u^/ (r s\n9) + fio an d the white arrows 
denote the meridional circulation. 

profiles have been obtained earlier from more specifically solar- 
like simulations (e.g., Brun & Toomre 2002; Brun et al. 2004; 
Brown et al. l2008l Ghizaru et al. [20T0t . 

In the slow rotation regime (Runs A1-A2) the kinetic en- 
ergy of meridional circulation and differential rotation are com- 
parable and comprise a few per cent of the total kinetic energy 
(columns 9 and 10 in Table [TJ. Increasing the Coriolis number 
further, increases the fraction of kinetic energy in the differen- 
tial rotation whereas that of the meridional circulation remains 
at first constant (Runs A3-A4), and finally drops close to zero 
(Runs A5-A6). In the three most rapidly rotating cases the dif- 
ferential rotation comprises more than 80 per cent of the total ki- 
netic energy. We also find that the meridional circulation shows a 
coherent pattern only for intermediate rotation rates (Runs A3- 
A5) where a single counter-clockwise cell per hemisphere ap- 
pears. In Run A6 the meridional flow is concentrated in a num- 
ber of small cells in accordance with earlier results (e.g., Miesch 
et al. 120001 Brun & Toomre 2002). We note that the rotation pro- 
files in Runs B3, CI, and D2 are similar to that in Run A6. 

The surface differential rotation of stars can be observa- 
tionally studied using photometric time series (e.g. Hall 1991) 
or with Doppler imaging methods (for a review, see Collier- 
Cameron 2007). The amount of surface differential rotation has 
been determined for some rapidly rotating pre- or main-sequence 
stars with varying spectral type (F, G, K, and M), systemati- 
cally showing solar-type differential rotation pattern with a faster 
equator and slower poles. The strength of the differential rotation 
shows a clear trend as function of the effective temperature, the 
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Fig. 17. Differential rotation parameter k^ according to Eq. (|46T > 
from Sets A (stars), B (diamonds), Run CI (cross), and Runs Dl 
and D2 (triangles). The dotted horizontal line indicates the zero 
level. 



shear being larger for hotter stars (see Fig. 1 of Collier-Cameron 
120071 1. Analysis of photometric time series, interpreting the pe- 
riod variations seen in the light curve analysis being due to dif- 
ferential rotation (e.g. Hall 119911 1. have established a relation 
An/n w fi-", with the values of n « 0.8-0.9. The observa- 
tional results are in rough agreement with theoretical predictions 
(e.g. Kitchatinov & Rudiger 1999), the theory predicting slightly 
weaker differential rotation in the rapid rotators than the actually 
observed values. 

We parameterise the differential rotation in our simulations 
with the quantity 



_ "eq 



pole 



a. 



An 



(46) 



where fl cq = fl(r 4 ,0 = 90°) and fi polo = 0(r 4 ,6> = 9 ). The 
results for the runs with Co ^ listed in Table Q] are shown in 
Fig. \T7\ We find that the anti-solar differential rotation peaks at 
Co « 1 and that ka turns positive for roughly Co ss 3. The val- 
ues in the rapid rotation (fcn ~ 0.2) end are comparable with the 
Sun (see also Chan 2010). It is not clear, however, how realistic 
it is to compare the current simulations with observations, i.e. 
even to argue that slowly rotating stars have anti-solar differen- 
tial rotation. It is clear that in the Sun the Coriolis number, and 
the radial length scale of convection, vary much more than in the 
current models so that it is not possible to reproduce equatorial 
acceleration and surface shear layer self-consistently in a sin- 
gle simulation. The situation may be different in slow rotators 
but observing their differential rotation is much more difficult. 
However, investigating the scaling of ka in the rapid rotation 
regime is likely worth pursuing (see also Brown et al. 2008). 



4. Conclusions 

The present results have demonstrated that the basic properties 
of Reynolds stress and turbulent heat flux found in Cartesian 
simulations are reproduced by simulations in spherical shells 
and wedges. This includes the signs of the off-diagonal compo- 
nents of Qij. In particular, the vertical stress, Q r $, is negative in 
both hemispheres when Co is small, but becomes positive near 
the top (and possibly also deeper down) when Co is large. This 
trend is well reproduced by the Cartesian simulations where Q xz 
is also negative for small Co, but becomes positive near the top 



when Co is large. These results coincide with earlier findings of 
Kapyla et al. (T2004] l. Chan d200Tb , and Robinson & Chan j20%T1 l. 

The horizontal stress Qe^, with the counterpart Q yz in the 
Cartesian model, is found to be positive in the northern hemi- 
sphere and have local maxima near the top and bottom of the 
domain. In spherical runs Qq^ is found to change sign near the 
poles for intermediate rotation. For rapid rotation, Q yz reaches 
a maximum near the top (or surface) around ±7° latitude - 
in agreement with earlier results (e.g., Chan 120011 Hupfer et 
al. 12005). We show that large-scale velocities due to the ba- 
nana cells near the equator are the main contribution to Q yz in 
Cartesian calculations. The spherical simulations reproduce such 
a sharp peak in the regime Co > 1, the peak being limited to a 
radially narrow region near the bottom of the domain. We find 
that the results for the Reynolds stress are weakly dependent on 
the Reynolds and Mach numbers. 

Furthermore, we find that Q r g is positive in the northern 
hemisphere, although for large values of Co the sign changes at 
the bottom of the convection zone. For the largest value of Co, 
Q t q is negative throughout the entire convection zone. A similar 
trend is seen in the Cartesian simulations, where Q xy is mostly 
positive but becomes negative near the bottom of the convection 
zone when rotation becomes strong enough, in accordance with 
Kapyla et al. J2004b 

The radial heat flux shows a strong dependence on latitude 
only when rotation is fairly rapid, i.e. Co > 1. This is associated 
with regions of the convection zone where the radial entropy gra- 
dient is decreased or even becomes positive. A partial explana- 
tion is that our setup (with a polytropic index of n — 1) is such 
that roughly 80 per cent of the energy is transported by radia- 
tive diffusion (cf. Brandenburg et al. 2005), making convection 
more easily suppressed than in a system where convection trans- 
ports a larger fraction. The latitudinal heat flux is equatorward 
for slow rotation and changes sign around Co « 1.A poleward 
heat flux is often used in breaking the Taylor-Proudman balance 
(e.g. Brandenburg et al. 1992). Longitudinal heat flux is mostly 
in the retrograde direction irrespective of the rotation rate. 

The turbulent heat conductivity \t is comparable to the first- 
order smoothing estimate with Strouhal number of the order of 
unity. The off-diagonal component xer is typically an order of 
magnitude smaller than the diagonal component \t in the rapid 
rotation regime. Similar results have been obtained previously 
from local convection simulations (e.g. Pulkkinen et al. [1993) 
and forced turbulence (Brandenburg et al. 2009). In mean-field 
models where anisotropic heat transport is invoked to break the 
Taylor-Proudman balance, the anisotropic part is typically of 
the same order of magnitude as the isotropic contribution (e.g. 
Brandenburg et al. 1992}. It is conceivable that the anisotropic 
contribution increases when the fraction of convective energy 
flux is increased. However, such a study is not within the scope 
of the present paper. 

As discussed in Sect. 13.11 the components of the Reynolds 
stress have contributions from diffusive and non-diffusive com- 
ponents. In future work we hope to be able to separate these two 
contributions, but in order to compare with earlier work, we have 
restricted ourselves to studying the components of the Reynolds 
stress directly. By making reasonable assumptions about the tur- 
bulent viscosity, it is indeed possible to obtain the relevant com- 
ponents of the A-effect, as was done by Pulkkinen et al. ( 119931 1. 
This is also true of global models, which also yield directly the 
global flow properties that can then be compared with corre- 
sponding mean field models, as was first done by Rieutord et 
al. ( [199 41). In a steady state, the Reynolds stress from the mean 
flow must then balance both the viscous stress and the Reynolds 
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stress from the fluctuations, as was demonstrated also by Miesch 
et al. (120081) . Such results are, however, dependent on the partic- 
ular properties of the model. 

In the present paper we find that in the slow and intermediate 
rotation regimes the differential rotation is anti-solar: the equator 
is rotating slower than the high latitudes. Such rotation profiles 
also coincide with the occurrence of coherent meridional circu- 
lation that is concentrated in a single counter-clockwise cell. In 
the rapid rotation regime, solar-like equatorial acceleration is ob- 
tained, but the differential rotation is confined to latitudes ±30° 
and the isocontours are aligned with the rotation vector. 

To reproduce the solar rotation profile at least two major ob- 
stacles remain. Firstly, the Taylor-Proudman balance must be 
broken. A possibility is to use subgrid-scale models where the 
present results for anisotropic heat transport can work as a guide. 
Secondly, the Coriolis number should decrease near the surface 
so that the transport of angular momentum is inward near the sur- 
face, leading to a surface shear layer as in the Sun. Here we can 
again introduce a subgrid-scale Reynolds stress guided by the 
present results. Studying such models, however, is postponed to 
future papers. 
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Appendix A: Dependence on domain size 

Above we have shown that we can recover many earlier results 
obtained in full spherical shells with wedges that span 150°in 
latitude and 90° in longitude. This gives at least a fourfold ad- 
vantage in terms of computation time in comparison to a full 
shell. However, it is important to study the range within which 
we can still recover the same results as with larger wedges. In 
order to study this we perform two additional sets of runs that 
are listed in Table I A. fl In Set E we vary the longitudinal extent 
from 22.5° to full 360°, with AO = 150° in all models. In Set F 
we keep the longitudinal extent fixed at Acf> = 90° and vary the 
latitudinal extent between 60° and 170°. As our base model we 
take Run A5 with fairly rapid rotation and complicated large- 
scale flows in the saturated state. 

Figure IaTI shows the latitudinal profiles of the off-diagonal 
components of the Reynolds stress from the middle of the con- 
vectively unstable layer and the rotation profiles as functions 
of radius from three latitudes from Set E and Run A5. The 
Reynolds stresses are very similar in the latitude range ±45° in 
runs with Aqb = 90° or larger. Somewhat larger differences are 
seen near the latitudinal boundaries. Runs El and E2 with the 
smallest longitude extents show the same qualitative behaviour 
for stress components Q r g and Qg^ but not for Q r< p. The rota- 
tion profiles for Runs A5, E3, and E4 with A(j> = 90° - 270° 
are very similar. The most obvious deviations from the trend oc- 
cur again for Runs El and E2 where the radial gradient of SI 
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Table A.l. Summary of the runs with varying A9 and A<fi. 



Run 


grid 


#i 


A8 


A(j> 


Ra 


Ma 


Re 


Co 


-Ether 


E 


kin 




EmcT / Ekiri 


-Grot /-Elan 




A5 


128 x 256 x 128 


15° 


150° 


90° 


3.1 


10° 


0.022 


36 


2.56 


0.111 


9.9 • 


10" 


r-r - 


0.002 


0.949 


-0.37 


El 


128 x 256 x 32 


15° 


150° 


22.5° 


3.1 


10 b 


0.019 


31 


2.91 


0.113 


7.0 ■ 


10" 




0.009 


0.963 


-0.43 


E2 


128 x 256 x 64 


15° 


150° 


45° 


3.1 


10 6 


0.020 


33 


2.77 


0.112 


7.7- 


10" 


-4 


0.008 


0.946 


-0.35 


E3 


128 x 256 x 256 


15° 


150° 


180° 


3.1 


10 6 


0.022 


37 


2.47 


0.113 


1.1 ■ 


10" 


-3 


0.002 


0.941 


-0.41 


E4 


128 x 256 x 384 


15° 


150° 


270° 


3.1 


10 6 


0.023 


39 


2.36 


0.111 


8.3 • 


10" 


-4 


0.002 


0.902 


-0.29 


E5 


128 x 256 x 512 


15° 


150° 


360° 


3.1 


10 6 


0.025 


41 


2.23 


0.112 


4.0 • 


10" 


-4 


0.001 


0.659 


-0.05 


Fl 


128 x 96 x 128 


60° 


60° 


90° 


3.1 


10° 


0.033 


31 


2.92 


0.115 


1.3 ■ 


10" 


-4 


0.001 


0.644 


+0.12 


F2 


128 x 160 x 128 


45° 


90° 


90° 


3.1 


10 6 


0.020 


32 


2.82 


0.114 


4.7- 


10" 


-4 


0.004 


0.899 


-0.13 


F3 


128 x 192 x 128 


30° 


120° 


90° 


3.1 


10 6 


0.021 


34 


2.68 


0.113 


4.3 • 


10" 


-4 


0.004 


0.885 


-0.08 


F4 


128 x 288 x 128 


5° 


170° 


90° 


3.1 


10 6 


0.020 


34 


2.71 


0.113 


3.7- 


10" 


-4 


0.006 


0.882 


-0.07 



is negative at the equator as opposed to the other runs where a 
positive gradient is found for r/R > 0.8. Surprisingly, Run E5 
with a full 360° longitude extent also deviates from the trend 
seen in the intermediate </>-extents: the qualitative trend of 
is similar but the magnitude of the differential rotation is re- 
duced. This is due to a non-axisymmetric m = 2 mode which 
is excited in this simulation. Large-scale hydrodynamical non- 
axisymmetries have been reported from rapidly rotating convec- 
tion (e.g. Brown et al. 2008 ). However, it is not clear whether the 
non-axisymmetry in our Run E5 is due to the same mechanism 
because of the slower rotation. 

Comparing simulations with different latitudinal extents 
(Fig. IA.2l l. we find that domains confined between ±45° latitude 
still reproduce the essential features of the solutions. This is par- 
ticularly clear for the Reynolds stresses which are very similar 
in the latitude range ±45° from the equator, with only Run Fl 
showing qualitatively different results in this range. There are 
also some differences at high latitudes between Runs A5 and 
F4. The rotation profiles are also very similar in the range ±30° 
with the exception of Run Fl. Run A5 also shows a deviating 
profile at high latitudes. 

These results suggest that a 90° longitude and 150° latitude 
extent is sufficient to capture the main features of the solutions 
at larger domains. The cost of this is that some features which 
are not of primary interest in the present study, such as the large- 
scale non-axisymmetric modes, are omitted. 
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Fig.A.l. Off-diagonal Reynolds stresses from the middle of the convection zone (upper row), and Q, as a function of radius at 
8 = 90° (lower row, left panel), 8 = 60° (middle panel), and 8 = 30° (right panel) for Runs E1-E5 and A5. Linestyles as indicated 
in the legend in the lower middle panel. 
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Fig.A.2. Same as Fig. lA.ll but for Runs F1-F4 and A5. The left panel on the lower row shows £1 from 8 = 45°. Linestyles as 
indicated in the legend in the lower left panel 
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